function [ y , t, out ] = simulaceDiff( model, linModel, inicializaceMPC, k_max, Ts, x0, u0, param )
%[ y , t, out ] = simulaceDiif( model, linModel, inicializaceMPC, k_max, Ts, x0, u0, param )
u=u0;
x(:,1)=x0;
outputSize=[];
for k=1:k_max;
    errx=x(:,k);
    save err errx u k
    [A,B,C,D]=linModel(x(:,k),u,param);
    sys=ss(A,B,C,D);
    if(isempty(outputSize))
        outputSize=size(C,1);
    end
    mpc=inicializaceMPC(c2d(sys,Ts),x(:,k),u,k, param);
    [du,up,vecp]=MPC(mpc);
    if(isempty(mpc.ref))
        out.ref(:,k)=zeros(outputSize,1);
    else
        out.ref(:,k)=mpc.ref(1:outputSize);
    end
    u=u+du;
    out.u(:,k)=u;
    out.up(:,:,k)=up;
    out.vecp(:,:,k)=vecp;
    [dState,y_k]=model(x(:,k),u,param);
    [dState2,~]=model(x(:,k)+Ts*dState,u,param);
    x(:,k+1)=x(:,k)+Ts/2*(dState+dState2);
    y(:,k)=y_k;
    t(k)=k*Ts;
end
out.x=x;
end

